## Preamble --------------------------------------------------------------------
library(arm)
library(multiwayvcov)
library(stargazer)
library(mlogit)
library(lfe)
library(ggplot2)
library(logitr)
library(tidyverse)
library(emdbook)
library(broom)
library(sf)
library(mapview)
library(tidycensus)
library(logitr)
library(ggmap) # you will need to use your own (free) API keys for both Google Maps and Stadia Maps and then use register_google() and register_stadiamaps() to make maps in this code file
library(tidytext)
library(SnowballC)

# load main data
wta_only <- read_rds("replication_data.rds") 
demog <- read_rds("demographics.rds")

## Summary statistics + descriptives tables ------------------------------------
(sensible <- nrow(subset(demog, housing_q3!= "" & !is.na(housing_q3)))) # 805
writeLines(as.character(round(sensible)), "Tables/n_sensible.tex",sep = "")

(sensible_wave_1 <- nrow(subset(demog, housing_q3!= "" & !is.na(housing_q3) & wave==1))) # 288
writeLines(as.character(round(sensible_wave_1)), "Tables/n_sensible_wave_1.tex",sep = "")

(sensible_wave_2 <- nrow(subset(demog, housing_q3!= "" & !is.na(housing_q3) & wave==2))) # 
writeLines(as.character(round(sensible_wave_2)), "Tables/n_sensible_wave_2.tex",sep = "")

(sensible_wave_3 <- nrow(subset(demog, housing_q3!= "" & !is.na(housing_q3) & wave==3))) # 
writeLines(as.character(round(sensible_wave_3)), "Tables/n_sensible_wave_3.tex",sep = "")

(demo_n <- nrow(subset(demog,  !is.na(demo_q10)))) # 589
writeLines(as.character(round(demo_n)), "Tables/n_demo.tex",sep = "")

completed <- demog %>%
  filter(!is.na(dev_prop_1_support)) %>%
  filter(!is.na(housing_q1)) # 748 when including individual proposals

(pct_ZIP <- mean(demog$Use.Zip.Method,na.rm=T))
writeLines(as.character(round(100*pct_ZIP)), "Tables/pct_ZIP.tex",sep = "")


descriptive_tab <- stargazer(as.data.frame(demog[,c("female",
                                      "white",
                                      "black",
                                      "latino",
                                      "asian",
                                      "age",
                                      "college",
                                      "income_gt90k",
                                      "own",
                                      "democrat",
                                      "liberal",
                                      "meeting")]),
                             summary = T, digits=2,
                             summary.stat=c("mean","sd","median","min","max","n"), 
                             covariate.labels = c("Female",
                                                  "White",
                                                  "Black",
                                                  "Latino",
                                                  "Asian",
                                                  "Age",
                                                  "College educated",
                                                  "Income \\textgreater 90k",
                                                  "Homeowner",
                                                  "Democrat",
                                                  "Liberal",
                                                  "Attended meeting"),
                             float = F,
                             out = "Tables/descriptive_tab_all.tex")

# wave 1
wave_1 <- subset(demog, wave==1)
descriptive_tab <- stargazer(as.data.frame(wave_1[,c("female",
                                       "white",
                                       "black",
                                       "latino",
                                       "asian",
                                       "age",
                                       "college",
                                       "income_gt90k",
                                       "own",
                                       "democrat",
                                       "liberal",
                                       "meeting")]), 
                             summary = T, digits=2,
                             summary.stat=c("mean","sd","median","min","max","n"), 
                             covariate.labels = c("Female",
                                                  "White",
                                                  "Black",
                                                  "Latino",
                                                  "Asian",
                                                  "Age",
                                                  "College educated",
                                                  "Income \\textgreater 90k",
                                                  "Homeowner",
                                                  "Democrat",
                                                  "Liberal",
                                                  "Attended meeting"),
                             float = F,
                             out = "Tables/descriptive_tab_wave1.tex")

# wave 2
wave_2 <- subset(demog, wave==2)
descriptive_tab <- stargazer(as.data.frame(wave_2[,c("female",
                                       "white",
                                       "black",
                                       "latino",
                                       "asian",
                                       "age",
                                       "college",
                                       "income_gt90k",
                                       "own",
                                       "democrat",
                                       "liberal",
                                       "meeting")]), 
                             summary = T, digits=2,
                             summary.stat=c("mean","sd","median","min","max","n"), 
                             covariate.labels = c("Female",
                                                  "White",
                                                  "Black",
                                                  "Latino",
                                                  "Asian",
                                                  "Age",
                                                  "College educated",
                                                  "Income \\textgreater 90k",
                                                  "Homeowner",
                                                  "Democrat",
                                                  "Liberal",
                                                  "Attended meeting"),
                             float = F,
                             out = "Tables/descriptive_tab_wave2.tex")

# wave 3
wave_3 <- subset(demog, wave==3)
descriptive_tab <- stargazer(as.data.frame(wave_3[,c("female",
                                       "white",
                                       "black",
                                       "latino",
                                       "asian",
                                       "age",
                                       "college",
                                       "income_gt90k",
                                       "own",
                                       "democrat",
                                       "liberal",
                                       "meeting")]), 
                             summary = T, digits=2,
                             summary.stat=c("mean","sd","median","min","max","n"), 
                             covariate.labels = c("Female",
                                                  "White",
                                                  "Black",
                                                  "Latino",
                                                  "Asian",
                                                  "Age",
                                                  "College educated",
                                                  "Income \\textgreater 90k",
                                                  "Homeowner",
                                                  "Democrat",
                                                  "Liberal",
                                                  "Attended meeting"),
                             float = F,
                             out = "Tables/descriptive_tab_wave3.tex")



## Results: support for proposals without compensation -------------------------
# broad cross-sectional support:
(pct_support <- nrow(subset(wta_only, support_bi == 1))/nrow(subset(wta_only, !is.na(support_bi))) )
writeLines(as.character(round(100*pct_support)), "Tables/pct_support.tex",sep = "")

# format clustering variables
wta_only$id <- as.factor(wta_only$id)
wta_only$proposal_id <- as.factor(wta_only$proposal_id)

# subset by homeownership status
wta_only_own <- subset(wta_only, own == 1)
wta_only_rent <- subset(wta_only, own == 0)



### Table D-6: support for proposals without compensation ----------------------
# generate models
lm1 <- felm(proposal_support ~ dist_km  + proposal_aff | 0 | 0 | id, wta_only)
lm2 <- felm(proposal_support ~ dist_km  + proposal_aff + own + income + white + black + latino + college + liberal + female + age + I(age^2) | 0 | 0 | id, wta_only)
lm3 <- felm(proposal_support ~ dist_km + proposal_aff*own + income + white + black + latino + college + liberal + female+ age + I(age^2) | 0 | 0 | id, wta_only)

stargazer(lm1, lm2, lm3,
          keep.stat = c("n", "rsq"),
          intercept.top = FALSE, intercept.bottom = TRUE,
          title = "Predictors of Support for Housing Proposals without Compensation",
          covariate.labels = c("Distance (km)", "Affordable", "Homeowner", "Income",
                               "White, non-Hispanic", "Black, non-Hispanic", "Hispanic", "College", "Liberal", "Female", 
                               "Age", "Age squared", "Affordable*Homeowner"),
          out = "Tables/ols_models.tex",
          label = "tab:ols_models",
          align = TRUE, type="text",
          dep.var.labels.include = F,
          float=F,
          dep.var.caption = "",
          star.cutoffs = c(.05, .01, .001),
          column.labels = c("\\shortstack{No \\\\ covariates}", "\\shortstack{With \\\\ covariates}",
                            "\\shortstack{Interact affordability \\\\ x homeownership}"),
          no.space = TRUE)

### Table H-10: Affordability interacted with ideology ---------------------------
# how does effect of affordability vary by ideology?
lm4a <- felm(proposal_support ~ dist_km + proposal_aff*liberal + own + income + white + black + latino + college + female + age + I(age^2) | 0 | 0 | id, wta_only)
lm4b <- felm(proposal_support ~ dist_km + proposal_aff*ideo + own + income + white + black + latino + college + female+ age + I(age^2) | 0 | 0 | id, wta_only)
lm4c <- felm(proposal_support ~ dist_km + proposal_aff*democrat + own + income + white + black + latino + college + female + age + I(age^2) | 0 | 0 | id, wta_only)
lm4d <- felm(proposal_support ~ dist_km + proposal_aff*pid + own + income + white + black + latino + college + female+ age + I(age^2) | 0 | 0 | id, wta_only)

stargazer(lm4a, lm4b, lm4c, lm4d,
          keep.stat = c("n", "rsq"),
          intercept.top = FALSE, intercept.bottom = TRUE,
          title = "Predictors of Support for Housing Proposals without Compensation",
          covariate.labels = c("Distance (km)", "Affordable", "Liberal", "Ideology", "Democrat", "Party ID",
                               "Homeowner", "Income",
                               "White, non-Hispanic", "Black, non-Hispanic", "Hispanic", "College", "Female",
                               "Age", "Age squared", "Affordable*Liberal", "Affordable*Ideology",
                               "Affordable*Democrat", "Affordable*Party ID"),
          out = "Tables/ols_models_aff_ideo.tex",
          label = "tab:ols_models_aff_ideo",
          align = TRUE, type="text",
          dep.var.labels.include = F,
          float=F,
          dep.var.caption = "",
          star.cutoffs = c(.05, .01, .001),
          no.space = TRUE)

### Table H-11: affordability, by survey wave ------------------------------------
# Analyses of support from affordability + covars, by wave
lm2a <- felm(proposal_support ~ dist_km  + proposal_aff + own + income + white + black + latino + college + liberal + female + age + I(age^2) | 0 | 0 | id, filter(wta_only,wave==1))
lm2b <- felm(proposal_support ~ dist_km  + proposal_aff + own + income + white + black + latino + college + liberal + female + age + I(age^2) | 0 | 0 | id, filter(wta_only,wave==2))
lm2c <- felm(proposal_support ~ dist_km  + proposal_aff + own + income + white + black + latino + college + liberal + female + age + I(age^2) | 0 | 0 | id, filter(wta_only,wave==3))

stargazer(lm2, lm2a, lm2b, lm2c,
          keep.stat = c("n", "rsq"),
          intercept.top = FALSE, intercept.bottom = TRUE,
          title = "Predictors of Support for Housing Proposals without Compensation",
          covariate.labels = c("Distance (km)", "Affordable", "Homeowner", "Income",
                               "White, non-Hispanic", "Black, non-Hispanic", "Hispanic", "College", "Liberal", "Female",
                               "Age", "Age squared"),
          out = "Tables/ols_models_wave.tex",
          label = "tab:ols_models_wave",
          align = TRUE, type="text",
          dep.var.labels.include = F,
          column.labels = c("All", "Wave 1", "Wave 2", "Wave 3"),
          dep.var.caption = "",
          float=F,
          star.cutoffs = c(.05, .01, .001),
          no.space = TRUE)



### Coefficient plots ---------
#### Figure 4: Effects on support for proposals without compensation ---------------
## without controls
coefs <- bind_rows(filter(tidy(lm1),term=="dist_km"),
                   filter(tidy(lm1),term=="proposal_aff")
) %>%
  mutate(term = factor(term,levels=c("proposal_aff","dist_km"),ordered=T))


## with controls
coefs_controls <- bind_rows(filter(tidy(lm2),term=="dist_km"),
                   filter(tidy(lm2),term=="proposal_aff")
) %>%
  mutate(term = factor(term,levels=c("proposal_aff","dist_km"),ordered=T))

coefs_both <- bind_rows(mutate(coefs_controls,controls="Yes"),
                        mutate(coefs,controls="No")
)

(coefplot_support_nocomp_both <- ggplot(coefs_both,aes(group=controls)) + 
    geom_errorbar(aes(x=term,ymin=(estimate + qnorm(0.025)*std.error), ymax=(estimate + qnorm(0.975)*std.error)), 
                  position=position_dodge(width=0.3),colour="black", width=0, size=0.5) +
    geom_errorbar(aes(x=term,ymin=(estimate + qnorm(0.05)*std.error), ymax=(estimate + qnorm(0.95)*std.error)), 
                  position=position_dodge(width=0.3),colour="black", width=0, size=1) +
    geom_point(aes(x=term,y=estimate,shape=controls),position=position_dodge(width=0.3),size=3) + 
    geom_hline(aes(yintercept=0),lty=2) + 
    scale_x_discrete("",breaks=c("dist_km","proposal_aff"),labels=c("Distance (in km)","Affordability")) + 
    scale_y_continuous("Effect on respondents' support") + 
    scale_shape_discrete("Demographic\ncontrols?",breaks=c("Yes","No")) + 
    coord_flip() + 
    theme_bw() + 
    theme(text = element_text(family="Arial Narrow"),legend.position = c(0.85,0.5),legend.text = element_text(lineheight=0.5),legend.spacing.y = unit(0.1,"cm"),legend.margin = margin(0,0,0,0))
)
ggsave(coefplot_support_nocomp_both,filename = paste0("Figures/coefplot_support_nocomp_both",".png"),height=3,width=4.5)


## Results: support for proposals with compensation ----------------------------
# average proposal support
writeLines(as.character(round(100*mean(wta_only$proposal_vote))), "Tables/pct_support_comp.tex",sep = "")

# full results
lm5 <- felm(proposal_vote ~ log(price) + dist_km + proposal_public + proposal_aff | 0| 0 | id , wta_only)
lm6 <- felm(proposal_vote ~ log(price) + dist_km + proposal_public + proposal_aff + own + income + white + black + latino + college + liberal  + female + age + I(age^2) | 0 | 0 | id , wta_only)
lm7 <- felm(proposal_vote ~ log(price)*proposal_aff + dist_km + proposal_public + own + income + white + black + latino + college + liberal  + female + age + I(age^2)| 0 | 0 | id , wta_only)
lm8 <- felm(proposal_vote ~ log(price)*proposal_public + dist_km + proposal_aff + own + income + white + black + latino + college + liberal  + female + age + I(age^2)| 0 | 0 | id , wta_only)


### Table D-7: support for proposals with compensation ----------------------
stargazer(lm5, lm6, lm7, lm8,
          keep.stat = c("n", "rsq"),
          intercept.top = FALSE, intercept.bottom = TRUE,
          title = "Predictors of Support for Housing Proposals with Compensation",
          covariate.labels = c("Compensation, logged", "Distance (km)", "Public benefits", "Affordable", "Homeowner",
                               "Income","White, non-Hispanic", "Black, non-Hispanic", "Hispanic", "College", 
                               "Democrat", "Female", "Age", "Age squared", "Compensation*Affordable", "Compensation*Public"),
          out = "Tables/ols_comp_models.tex",
          label = "tab:ols_comp_models",
          align = TRUE, type="text",
          float=F,
          dep.var.labels.include = F,
          dep.var.caption = "",
          star.cutoffs = c(.05, .01, .001),
          column.labels = c("\\shortstack{No \\\\ covariates}", "\\shortstack{With \\\\ covariates}",
                            "\\shortstack{Interact price \\\\ x affordability}", "\\shortstack{Interact price \\\\ x form of comp.}"),
          no.space = TRUE)

### Table H-12: effects of compensation by survey wave -------------------------

# effect of compensation by wave
lm6a <- felm(proposal_vote ~ log(price) + dist_km + proposal_public + proposal_aff + own + income + white + black + latino + college + liberal  + female + age + I(age^2) | 0 | 0 | id , filter(wta_only,wave==1))
lm6b <- felm(proposal_vote ~ log(price) + dist_km + proposal_public + proposal_aff + own + income + white + black + latino + college + liberal  + female + age + I(age^2) | 0 | 0 | id , filter(wta_only,wave==2))
lm6c <- felm(proposal_vote ~ log(price) + dist_km + proposal_public + proposal_aff + own + income + white + black + latino + college + liberal  + female + age + I(age^2) | 0 | 0 | id , filter(wta_only,wave==3))

stargazer(lm6, lm6a, lm6b, lm6c,
          keep.stat = c("n", "rsq"),
          intercept.top = FALSE, intercept.bottom = TRUE,
          title = "Predictors of Support for Housing Proposals with Compensation",
          covariate.labels = c("Compensation, logged", "Distance (km)", "Public benefits", "Affordable", "Homeowner",
                               "Income","White, non-Hispanic", "Black, non-Hispanic", "Hispanic", "College",
                               "Democrat", "Female", "Age", "Age squared"),
          out = "Tables/ols_comp_models_wave.tex",
          label = "tab:ols_comp_models_wave",
          align = TRUE, type="text",
          dep.var.labels.include = F,
          float=F,
          dep.var.caption = "",
          star.cutoffs = c(.05, .01, .001),
          column.labels = c("All", "Wave 1",
                            "Wave 2", "Wave 3"),
          no.space = TRUE)



### Table D-8: effects of compensation by homeowner vs. renter -------------------------
comp_support_own <- felm(proposal_vote ~ log(price)*proposal_aff + proposal_public| 0 | 0 | id, wta_only_own)
comp_support_own_c <- felm(proposal_vote ~ log(price)*proposal_aff + proposal_public + dist_km + income + white + black + latino + college + liberal  + female  + age + I(age^2)| 0 | 0 | id, wta_only_own)

comp_support_rent <- felm(proposal_vote ~ log(price)*proposal_aff + proposal_public| 0 | 0 | id, wta_only_rent)
comp_support_rent_c <- felm(proposal_vote ~ log(price)*proposal_aff + proposal_public + dist_km + income + white + black + latino + college + liberal  + female  + age + I(age^2) | 0 | 0 | id, wta_only_rent)


stargazer(comp_support_rent, comp_support_rent_c, comp_support_own, comp_support_own_c,
          keep.stat = c("n", "rsq"),
          intercept.top = FALSE, intercept.bottom = TRUE,
          title = "Predictors of Support for Housing Proposals with Compensation",
          covariate.labels = c("Compensation, logged", "Affordable", "Public benefits" ,"Distance (km)", "Income",
                               "White, non-Hispanic", "Black, non-Hispanic",  "Hispanic", "College",
                               "Liberal", "Female", "Age", "Age squared", "Compensation*Affordable", "Constant"),
          out = "Tables/ols_models_comp_own.tex",
          label = "tab:ols_models_comp_own",
          align = TRUE, type="text",
          dep.var.labels.include = F,
          float = F,
          dep.var.caption = "",
          star.cutoffs = c(.05, .01, .001),
          column.labels = c("\\shortstack{Renters: \\\\ No \\\\ covariates}", "\\shortstack{Renters: \\\\ With \\\\ covariates}",
                            "\\shortstack{Homeowners: \\\\ No \\\\ covariates}", "\\shortstack{Homeowners: \\\\ With \\\\ covariates}"),
          no.space = TRUE)


### Table D-9: effects of compensation by income -------------------------------
# Does effect of compensation vary by income?
comp_support_own_inc <- felm(proposal_vote ~ log(price)*income + proposal_public| 0 | 0 | id, subset(wta_only_own, proposal_aff == 0))
comp_support_own_inc_c <- felm(proposal_vote ~ log(price)*income + proposal_public + dist_km  + white + black + latino + college + liberal  + female  + age + I(age^2)| 0 | 0 | id, subset(wta_only_own, proposal_aff == 0))

comp_support_rent_inc <- felm(proposal_vote ~ log(price)*income + proposal_public| 0 | 0 | id, subset(wta_only_rent, proposal_aff == 0))
comp_support_rent_inc_c <- felm(proposal_vote ~ log(price)*income + proposal_public + dist_km + white + black + latino + college + liberal  + female  + age + I(age^2) | 0 | 0 | id, subset(wta_only_rent, proposal_aff == 0))

stargazer(comp_support_rent_inc, comp_support_rent_inc_c, comp_support_own_inc, comp_support_own_inc_c,
          keep.stat = c("n", "rsq"),
          intercept.top = FALSE, intercept.bottom = TRUE,
          title = "Predictors of Support for Housing Proposals with Compensation",
          covariate.labels = c("Compensation, logged", "Income",  "Public benefits", "Distance (km)", 
                               "White, non-Hispanic", "Black, non-Hispanic",  "Hispanic", "College", 
                               "Liberal", "Female", "Age", "Age squared", "Compensation*Income"),
          out = "Tables/ols_models_comp_own_inc.tex",
          label = "tab:ols_models_comp_own_inc",
          align = TRUE, type="text",
          dep.var.labels.include = F,
          float=F,
          dep.var.caption = "",
          star.cutoffs = c(.05, .01, .001),
          column.labels = c("\\shortstack{Renters: \\\\ No \\\\ covariates}", "\\shortstack{Renters: \\\\ With \\\\ covariates}",
                            "\\shortstack{Homeowners: \\\\ No \\\\ covariates}", "\\shortstack{Homeowners: \\\\ With \\\\ covariates}"),
          no.space = TRUE)



### Nonparametric visuals ------------------------------------------------------

#### Figure D-1a and D-1b ------------------------------------------------------

## Renters only

# affordable
png("Figures/comp_support_rent_aff.png", width=320, height = 320, units = "px", pointsize = 14)
ggplot(subset(wta_only_rent, proposal_aff == 1), aes(x=log(price), y=proposal_vote))+ 
  geom_jitter(size=2, alpha=0.2, position = position_jitter(height = .02))+
  geom_smooth(method="loess", colour="blue", size=1.5)+
  xlab("Compensation ($, logged)")+
  ylab("Probability of Support")+
  theme_bw() 
dev.off()

# market rate
png("Figures/comp_support_rent_mar.png", width=320, height = 320, units = "px", pointsize = 14)
ggplot(subset(wta_only_rent, proposal_aff == 0), aes(x=log(price), y=proposal_vote))+ geom_jitter(size=2, alpha=0.2, position = position_jitter(height = .02))+
  geom_smooth(method="loess", colour="blue", size=1.5)+
  xlab("Compensation ($, logged)")+
  ylab("Probability of Support")+
  theme_bw() 
dev.off()

#### Figure D-2a and D-2b ------------------------------------------------------

## Homeowners only

# affordable
png("Figures/comp_support_own_aff.png", width=320, height = 320, units = "px", pointsize = 14)
ggplot(subset(wta_only_own, proposal_aff == 1), aes(x=log(price), y=proposal_vote))+ 
  geom_jitter(size=2, alpha=0.2, position = position_jitter(height = .02))+
  geom_smooth(method="loess", colour="blue", size=1.5)+
  xlab("Compensation ($, logged)")+
  ylab("Probability of Support")+
  theme_bw()
dev.off()

# market rate
png("Figures/comp_support_own_mar.png", width=320, height = 320, units = "px", pointsize = 14)
ggplot(subset(wta_only_own, proposal_aff == 0), aes(x=log(price), y=proposal_vote))+ 
  geom_jitter(size=2, alpha=0.2, position = position_jitter(height = .02))+
  geom_smooth(method="loess", colour="blue", size=1.5)+
  xlab("Compensation ($, logged)")+
  ylab("Probability of Support")+
  theme_bw()
dev.off()


### Coefficient plots ------------------------------------------------------
#### Figure 5: Effects on support for proposals with compensation ---------------
## no controls
coefs <- bind_rows(filter(tidy(lm5),term=="log(price)"),
                   filter(tidy(lm5),term=="dist_km"),
                   filter(tidy(lm5),term=="proposal_aff"),
                   filter(tidy(lm5),term=="proposal_public")
) %>%
  mutate(term = factor(term,levels=c("log(price)","proposal_public","proposal_aff","dist_km"),ordered=T))

## with controls
coefs_controls <- bind_rows(filter(tidy(lm6),term=="log(price)"),
                            filter(tidy(lm6),term=="dist_km"),
                            filter(tidy(lm6),term=="proposal_aff"),
                            filter(tidy(lm6),term=="proposal_public")
) %>%
  mutate(term = factor(term,levels=c("log(price)","proposal_public","proposal_aff","dist_km"),ordered=T))

coefs_both <- bind_rows(mutate(coefs_controls,controls="Yes"),
                        mutate(coefs,controls="No")
)

(coefplot_support_comp_both <- ggplot(coefs_both,aes(group=controls)) + 
    geom_errorbar(aes(x=term,ymin=(estimate + qnorm(0.025)*std.error), ymax=(estimate + qnorm(0.975)*std.error)), 
                  position=position_dodge(width=0.3),colour="black", width=0, size=0.5) +
    geom_errorbar(aes(x=term,ymin=(estimate + qnorm(0.05)*std.error), ymax=(estimate + qnorm(0.95)*std.error)), 
                  position=position_dodge(width=0.3),colour="black", width=0, size=1) +
    geom_point(aes(x=term,y=estimate,shape=controls),position=position_dodge(width=0.3),size=3) + 
    geom_hline(aes(yintercept=0),lty=2) + 
    scale_x_discrete("",breaks=c("log(price)","proposal_public","dist_km","proposal_aff"),labels=c("log(Price in $)","Public vs. private benefits","Distance (in km)","Affordability")) + 
    scale_y_continuous("Effect on respondents' binary support") + 
    scale_shape_discrete("Demographic\ncontrols?",breaks=c("Yes","No")) + 
    coord_flip() + 
    theme_bw() + 
    theme(text = element_text(family="Arial Narrow"),legend.position = c(0.85,0.18),legend.text = element_text(lineheight=0.5),legend.spacing.y = unit(0.1,"cm"),legend.margin = margin(0,0,0,0))
)
ggsave(coefplot_support_comp_both,filename = paste0("Figures/coefplot_support_comp_both",".png"),height=4,width=4.5)


#### Figure 6: Effects on support for proposals with compensation, interacted with affordability -----------

## with interaction:
coefs_controls_intx <- tibble(
  term=c("dist_km","proposal_aff","proposal_public","log(price)","log(price)"),
  subgroup=c("All Proposals","All Proposals","All Proposals","Market-Rate","Affordable"),
  estimate=c(
             filter(tidy(lm7),term=="dist_km")$estimate,
             filter(tidy(lm7),term=="proposal_aff")$estimate,
             filter(tidy(lm7),term=="proposal_public")$estimate,
             filter(tidy(lm7),term=="log(price)")$estimate,
             filter(tidy(lm7),term=="log(price):proposal_aff")$estimate+filter(tidy(lm7),term=="log(price)")$estimate),
  std.error=c(filter(tidy(lm7),term=="dist_km")$std.error,
              filter(tidy(lm7),term=="proposal_aff")$std.error,
              filter(tidy(lm7),term=="proposal_public")$std.error,
              filter(tidy(lm7),term=="log(price)")$std.error,
              sqrt(filter(tidy(lm7),term=="log(price):proposal_aff")$std.error^2+filter(tidy(lm7),term=="log(price)")$std.error^2))
) %>%
  mutate(term = factor(term,levels=c("log(price)","proposal_public","proposal_aff","dist_km"),ordered=T))



(coefplot_support_comp_intx <- ggplot(coefs_controls_intx,aes(group=subgroup)) + 
    geom_errorbar(aes(x=term,ymin=(estimate + qnorm(0.025)*std.error), ymax=(estimate + qnorm(0.975)*std.error)), 
                  position=position_dodge(width=0.3),colour="black", width=0, size=0.5) +
    geom_errorbar(aes(x=term,ymin=(estimate + qnorm(0.05)*std.error), ymax=(estimate + qnorm(0.95)*std.error)), 
                  position=position_dodge(width=0.3),colour="black", width=0, size=1) +
    geom_point(aes(x=term,y=estimate,shape=subgroup),position=position_dodge(width=0.3), size=3) + 
    geom_hline(aes(yintercept=0),lty=2) + 
    annotate(geom="text",x=1.5,y=0.2,label="Market-Rate",family="Arial Narrow",hjust=0) + 
    geom_curve(aes(x = 1.5, y = 0.19, xend = 1.2, yend = 0.06),arrow = arrow(length = unit(0.01, "npc"),type="closed"), size = 0.4,curvature = 0.3,angle=90) + 
    annotate(geom="text",x=0.5,y=0.2,label="Affordable",family="Arial Narrow",hjust=0) + 
    geom_curve(aes(x = 0.5, y = 0.19, xend = 0.8, yend = 0.025),arrow = arrow(length = unit(0.01, "npc"),type="closed"), size = 0.4,curvature = -0.3,angle=90) + 
    scale_x_discrete("",breaks=c("log(price)","proposal_public","dist_km","proposal_aff"),labels=c("log(Price in $)","Public vs. private benefits","Distance (in km)","Affordability")) + 
    scale_y_continuous("Effect on respondents' binary support") + 
    scale_shape_discrete("Proposal Type",breaks=c("All Proposals","Market-Rate","Affordable")) +
    coord_flip() + 
    theme_bw() + 
    theme(text = element_text(family="Arial Narrow"),legend.position = "none")
)
ggsave(coefplot_support_comp_intx,filename = paste0("Figures/coefplot_support_comp_intx",".png"),height=4,width=5.5)




# Appendix E: Analyses using Multinomial & Mixed Logit -------------------------
wta_aff <- read_rds("replication_data_logitr.rds")

# remove any missing data and add interaction variables for logitr
logitr_data <- wta_aff %>%
  filter(!is.na(own)) %>%
  filter(!is.na(proposal_aff)) %>%
  filter(!is.na(build)) %>%
  filter(!is.na(choiceID)) %>%
  filter(!is.na(price)) %>%
  mutate(build_own = build*own, 
         price_own = price*own, 
         aff_own = proposal_aff*own, 
         price_aff = price*proposal_aff, 
         build_aff = build*proposal_aff,
         obsID = choiceID, 
         choice = proposal_vote, 
         price_dist = price*dist_2sd, 
         build_dist = build*dist_2sd) %>%
  arrange(obsID)

## H5a Compensation increases support
# index as number
logitr_data$obs_index <- rep(1:(nrow(logitr_data)/2), each = 2)
length(unique(logitr_data$id)) # 397 unique respondents here


### Effect of compensation by proximity of proposal ----------------------------
## mnl
summary(logitr_data$dist_km)

# near
logitr_data_near <- subset(logitr_data, dist_km < .54)
mnl_h1_near <- logitr(
  data = logitr_data_near,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)
summary(mnl_h1_near)

logit1 <- glm(formula = choice ~ price*build,
              data = logitr_data_near,
              family = binomial)
summary(logit1)$coefficients


# far
logitr_data_far <- subset(logitr_data, dist_km > .54)
mnl_h1_far <- logitr(
  data = logitr_data_far,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)

new_data <- data.frame("obs_index" = c(1,1,2,2,3,3,4,4,5,5), 
                       "price" = rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1), 
                       "build" = rep(c(0, 1), 5))

probs_near <- predict(
  mnl_h1_near,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_near$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_near$Proximity <- rep("Near", 10)
probs_near <- subset(probs_near, price > 0)

probs_far <- predict(
  mnl_h1_far,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_far$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_far$Proximity <- rep("Far", 10)
probs_far <- subset(probs_far, price > 0)

probs <- rbind(probs_near, probs_far)


#### Figure E3 ---------
(logit_h1 <- ggplot(probs, aes(x = price, y = predicted_prob, color = Proximity)) + 
  geom_line(stat = "identity") + 
  geom_ribbon(aes(ymin = predicted_prob_lower, ymax = predicted_prob_upper), alpha = 0.1) +
  scale_y_continuous(limits = c(0, 1)) + labs(x = "Compensation", y = "Expected Support Probability")+
  theme_bw()
)
ggsave(logit_h1,filename = paste0("Figures/logit_h1",".png"),height=5,width=5)


### Effect of compensation by affordability of proposal ------------------------
### h2 price * aff
`%!in%` = Negate(`%in%`)
table(logitr_data$proposal_aff)
aff <- subset(logitr_data, proposal_aff == 1)
id_aff <- subset(logitr_data, obs_index %in% aff$obs_index)
id_mar <- subset(logitr_data, obs_index %!in% aff$obs_index)

# market rate
mnl_h1_mar <- logitr(
  data = id_mar,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)
summary(mnl_h1_mar)

# affordable
mnl_h1_aff <- logitr(
  data = id_aff,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)
summary(mnl_h1_aff)

new_data <- data.frame("obs_index" = c(1,1,2,2,3,3,4,4,5,5), 
                       "price" = rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1), 
                       "build" = rep(c(0, 1), 5))

probs_mar <- predict(
  mnl_h1_mar,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_mar$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_mar$Type <- rep("Market-rate", 10)
probs_mar <- subset(probs_mar, price > 0)

probs_aff <- predict(
  mnl_h1_aff,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_aff$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_aff$Type <- rep("Affordable", 10)
probs_aff <- subset(probs_aff, price > 0)

probs <- rbind(probs_mar, probs_aff)


#### Figure E-4 ------------
(logit_h2 <- ggplot(probs, aes(x = price, y = predicted_prob, color = Type)) + 
  geom_line(stat = "identity") + 
  geom_ribbon(aes(ymin = predicted_prob_lower, ymax = predicted_prob_upper), alpha = 0.1) +
  scale_y_continuous(limits = c(0, 1)) + labs(x = "Compensation", y = "Expected Support Probability")+
  theme_bw()
)
ggsave(logit_h2, filename = paste0("Figures/logit_h2",".png"), height=5, width=5)


### Effect of compensation by homeownership ------------------------------------
# h3 Effect of price by homeownership level
owners <- subset(logitr_data, own == 1)
id_own <- subset(logitr_data, obs_index %in% owners$obs_index)
id_rent <- subset(logitr_data, obs_index %!in% owners$obs_index)

# homeowners
mnl_h1_own <- logitr(
  data = id_own,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)
summary(mnl_h1_own)

# renters
mnl_h1_rent <- logitr(
  data = id_rent,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)
summary(mnl_h1_rent)

test <- logitr(
  data = logitr_data,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build", "price_own", "build_own")
)
summary(test)


new_data <- data.frame("obs_index" = c(1,1,2,2,3,3,4,4,5,5), 
                       "price" = rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1), 
                       "build" = rep(c(0, 1), 5))
new_data

probs_own <- predict(
  mnl_h1_own,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_own$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_own$Tenure <- rep("Homeowners", 10)
probs_own <- subset(probs_own, price > 0)

probs_rent <- predict(
  mnl_h1_rent,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_rent$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_rent$Tenure <- rep("Renters", 10)
probs_rent <- subset(probs_rent, price > 0)

probs <- rbind(probs_rent, probs_own)
probs
probs$Tenure <- factor(probs$Tenure, levels = c("Renters", "Homeowners"))


#### Figure E-5 -----------------------------------------------------------------
logit_h3 <- ggplot(probs, aes(x = price, y = predicted_prob, color = Tenure)) + 
  geom_line(stat = "identity") + 
  geom_ribbon(aes(ymin = predicted_prob_lower, ymax = predicted_prob_upper), alpha = 0.1) +
  scale_y_continuous(limits = c(0, 1)) + labs(x = "Compensation", y = "Expected Support Probability")+
  theme_bw()
logit_h3
ggsave(logit_h3, filename = paste0("Figures/logit_h3",".png"), height=5, width=5)


### Effect of compensation by affordability of proposals, by homeownership of respondent ---------
# h4 renters only for affordable
aff <- subset(logitr_data, proposal_aff == 1)
id_aff <- subset(logitr_data, own == 0 & obs_index %in% aff$obs_index)
id_mar <- subset(logitr_data, own == 0 & obs_index %!in% aff$obs_index)

# market rate
mnl_h1_mar <- logitr(
  data = id_mar,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)
summary(mnl_h1_mar)

# affordable
mnl_h1_aff <- logitr(
  data = id_aff,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)
summary(mnl_h1_aff)

new_data <- data.frame("obs_index" = c(1,1,2,2,3,3,4,4,5,5), 
                       "price" = rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1), 
                       "build" = rep(c(0, 1), 5))
new_data

probs_mar <- predict(
  mnl_h1_mar,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_mar$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_mar$Affordability <- rep("Market-rate", 10)
probs_mar <- subset(probs_mar, price > 0)

probs_aff <- predict(
  mnl_h1_aff,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_aff$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_aff$Affordability <- rep("Affordable", 10)
probs_aff <- subset(probs_aff, price > 0)

probs <- rbind(probs_mar, probs_aff)
probs


#### Figure E-6: among renters -------------------------------------------------
logit_h4 <- ggplot(probs, aes(x = price, y = predicted_prob, color = Affordability)) + 
  geom_line(stat = "identity") + 
  geom_ribbon(aes(ymin = predicted_prob_lower, ymax = predicted_prob_upper), alpha = 0.1) +
  scale_y_continuous(limits = c(0, 1)) + labs(x = "Compensation", y = "Expected Support Probability")+
  theme_bw()
logit_h4
ggsave(logit_h4, filename = paste0("Figures/logit_h4",".png"), height=5, width=5)


### Effect of compensation by form of compensation ----------------------------

# h5 effect of form of compensation
public <- subset(logitr_data, proposal_public == 1)
id_pub <- subset(logitr_data, obs_index %in% public$obs_index)
id_priv <- subset(logitr_data, obs_index %!in% public$obs_index)

# public
mnl_h1_pub <- logitr(
  data = id_pub,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)
summary(mnl_h1_pub)

# private
mnl_h1_priv <- logitr(
  data = id_priv,
  outcome = "choice",
  obsID = "obs_index",
  pars = c("price", "build")
)
summary(mnl_h1_priv)

new_data <- data.frame("obs_index" = c(1,1,2,2,3,3,4,4,5,5), 
                       "price" = rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1), 
                       "build" = rep(c(0, 1), 5))
new_data

probs_pub <- predict(
  mnl_h1_pub,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_pub$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_pub$Form <- rep("Public", 10)
probs_pub <- subset(probs_pub, price > 0)

probs_priv <- predict(
  mnl_h1_priv,
  newdata = new_data,
  obsID   = "obs_index",
  ci      = 0.95, returnData = TRUE
)
probs_priv$price <- rep(c(0, 100, 0, 2500, 0, 5000, 0, 7500, 0, 10000), 1)
probs_priv$Form <- rep("Private", 10)
probs_priv <- subset(probs_priv, price > 0)

probs <- rbind(probs_pub, probs_priv)
probs

#### Figure E-7 ---------------
(logit_h5 <- ggplot(probs, aes(x = price, y = predicted_prob, color = Form)) + 
  geom_line(stat = "identity") + 
  geom_ribbon(aes(ymin = predicted_prob_lower, ymax = predicted_prob_upper), alpha = 0.1) +
  scale_y_continuous(limits = c(0, 1)) + labs(x = "Compensation", y = "Expected Support Probability")+
  theme_bw()
)
ggsave(logit_h5, filename = paste0("Figures/logit_h5",".png"), height=5, width=5)



## Analysis of open-ended comments ---------------------------------------------
df <- read_rds("open_ended_text.rds")

# market-rate
market <- subset(df, dev_prop_1_afford_treat == "The units would be rented at whatever price the local market supports."  )
market <- tibble(id = market$user.id, text = market$dev_prop_1_cond_comments)
market <- subset(market, !is.na(text))
market <- market %>%
  unnest_tokens(word, text) %>%
  mutate(stem = wordStem(word)) %>%
  anti_join(stop_words) %>% # remove stop words
  filter(!grepl("[[:digit:]]", stem)) # remove numbers

market <- market %>%
  count(stem, sort = TRUE) %>%
  mutate(total = sum(n)) %>%
  mutate(tf = n/total) %>%
  arrange(desc(n))

market <- market %>%
  rename(tf_market = tf, n_market = n, total_market = total)

# affordable
affordable <- subset(df, dev_prop_1_afford_treat == "Half of the units would be occupied by low-income housing voucher recipients." )
affordable <- tibble(id = affordable$user.id, comments = affordable$dev_prop_1_cond_comments)
affordable <- subset(affordable, !is.na(comments))
affordable <- affordable %>%
  unnest_tokens(word, comments) %>%
  mutate(stem = wordStem(word)) %>%
  anti_join(stop_words) %>%
  filter(!grepl("[[:digit:]]", stem))

affordable <- affordable %>%
  count(stem, sort = TRUE) %>%
  mutate(total = sum(n)) %>%
  mutate(tf = n/total) %>%
  arrange(desc(n))

affordable <- affordable %>%
  rename(tf_affordable = tf, n_affordable = n, total_affordable = total)

stem_tf <- left_join(affordable, market, by = c("stem" = "stem")) %>%
  mutate(ratio = (tf_affordable/tf_market)) %>%
  mutate(log2_ratio = log2(tf_affordable/tf_market))

# select top stems
stem_tf_affordable_top <- stem_tf %>%
  top_n(n=17, n_affordable) 

stem_tf_market_top <- stem_tf %>%
  top_n(n=17, n_market)

stem2 <- rbind(stem_tf_affordable_top, stem_tf_market_top)
dupes <- which(duplicated(stem2$stem))

# remove duplicate entry to get slim list
if(length(dupes)>0) {
  stem2 <- stem2[-dupes, ]
}

# 21 total terms now, due to duplicates
stem2_affordable <- stem2 %>%
  top_n(n=7, log2_ratio)

stem2_market <- stem2 %>%
  top_n(n=14, -log2_ratio)

stem3 <- rbind(stem2_affordable, stem2_market)
stem_tf_small <- stem3 %>%
  distinct()


(p <- stem_tf_small %>%
    mutate(stem = reorder(stem, log2_ratio)) %>%
    arrange(desc(log2_ratio)) %>%
    ggplot(aes(stem, log2_ratio) ) +
    ylab(label = "Log2 ratio of term frequencies in comments\n(Affordable treatment / Market-rate treatment)") +
    xlab(label = "") +
    geom_col(show.legend = FALSE) +
    coord_flip() + 
    theme_bw() 
)

ggsave(p,filename = "Figures/tf_ratios_affordability.png", width=6, height=5)

## Maps of proposals/respondent locations -----------
a_sf <- read_rds("respondent_locations.rds")

props <- read_rds("proposal_locations.rds")

# mapview(a_sf, zcol = "code", legend = F, zoom = 12)
# ggsave("Figures/respondents.pdf", height = 4, length = 4 )

bos_map <- ggmap::get_map(location=c(lon = -71.075, lat = 42.324), zoom = 12,source = "stadia",maptype="alidade_smooth")

### Figure 3a -------------------------------------------------------------------
(respondent_map <- ggmap(bos_map, extent = "device") + 
   geom_sf(data = a_sf,  color="#414487FF",alpha=0.6, inherit.aes = F)  +
   scale_color_viridis_c() + 
   theme_bw() + 
   theme(legend.justification=c(0,0), legend.position=c(0.7,0.025),
         axis.line = element_blank(), # get rid of axis lines and long/lat markers
         axis.text = element_blank(),
         axis.ticks = element_blank(),
         axis.title = element_blank())
)
ggsave("Figures/respondent_map.png",respondent_map, height = 4, width = 4)

### Figure 3b -------------------------------------------------------------------
(proposal_map <- ggmap(bos_map, extent = "device") + 
   geom_sf(data = props,  color="red",alpha=0.6, inherit.aes = F)  +
   scale_color_viridis_c() + 
   theme_bw() + 
   theme(legend.justification=c(0,0), legend.position=c(0.7,0.025),
         axis.line = element_blank(), # get rid of axis lines and long/lat markers
         axis.text = element_blank(),
         axis.ticks = element_blank(),
         axis.title = element_blank())
)
ggsave("Figures/proposal_map.png", proposal_map, height = 4, width = 4)




## Boston CBAs -------
df3 <- read_rds("boston_cbas.rds")

boston <- read_rds("boston_tracts.rds")

bos_map <- ggmap::get_map(location=c(lon = -71.075, lat = 42.324), zoom = 12)

### Figure 1 -------------------------------------------------------------------
(income_cba <- ggmap(bos_map, extent = "device") + 
    geom_sf(subset(boston, !is.na(med_incomeE) & popE>22), # removed low-pop tracts to get rid of mostly-park tracts
            mapping = aes(fill = med_incomeE), lwd=0,inherit.aes = F) +
    scale_fill_gradient(high = "#e34a33", low = "#fee8c8", guide = "colorbar",labels=scales::dollar_format()) +
    theme_bw() + 
    guides(fill=guide_colorbar(title="Household Median Income")) + 
    geom_sf(data = df3,  color = "black", inherit.aes = F)  +
    theme(legend.justification=c(0,0), legend.position=c(0.7,0.025),
          axis.line = element_blank(), # get rid of axis lines and long/lat markers
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          axis.title = element_blank())
)

ggsave("Figures/income_cba.pdf",income_cba, height = 8, width = 8)
